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This paper presents a Kalman filter using a seven-component attitude state vector 
comprising the angular momentum components in an inertial reference frame, the angular 
momentum components in the body frame, and a rotation angle. The relatively slow 
variation of these parameters makes this parameterization advantageous for spinning 
spacecraft attitude estimation. The filter accounts for the constraint that the magnitude of 
the angular momentum vector is the same in the inertial and body frames by employing a 
reduced six-component error state. Four variants of the filter, defined by different choices 
for the reduced error state, are tested against a quaternion-based filter using simulated data 
for the THEMIS mission. Three of these variants choose three of the components of the 
error state to be the infinitesimal attitude error angles, facilitating the computation of 
measurement sensitivity matrices and causing the usual 3x3 attitude covariance matrix to be 
a submatrix of the 6x6 covariance of the error state. These variants differ in their choice for 
the other three components of the error state. The variant employing the infinitesimal 
attitude error angles and the angular momentum components in an inertial reference frame 
as the error state shows the best combination of robustness and efficiency in the simulations. 

Attitude estimation results using THEMIS flight data are also presented. 

I. Introduction 

Attitude estimation is often more difficult for spinning spacecraft than for three-axis stabilized spacecraft. The 
parameters representing the spacecraft attitude and its time rate of change vary more rapidly in the spinning case, 
and gyro rate measurements are often lacking, requiring Euler’s equations for modeling the attitude dynamics. This 
paper uses a seven-parameter angular-momentum-based representation that is advantageous for this application. 1 
The seven state vector elements are the angular momentum components in an inertial reference frame, the angular 
momentum components in the spacecraft’s body frame, and a rotation angle. These parameters are subject to the 
constraint that the magnitude of the angular momentum vector is the same in the inertial and body frames. 

We have developed a series of extended Kalman filters (EKFs) employing this representation, to which we give the 
generic name SpinKF. The constraint on the state vector allows us to employ a six-component error state instead of 
the error vector of the full seven-component state, in parallel with the procedure commonly used to estimate the 
constrained four-component quaternion representation of attitude. 2 The conceptual advantages of this dimensional 
reduction, as more truly representing the actual degrees of freedom of the system, have been debated at length, 3-9 but 
the computational advantages are indisputable. The different versions of SpinKF all use the same seven-component 
state but differ in their specification of the six-component error state. The first, SpinKFl, 10 was incorporated into the 
attitude ground support system for spinning spacecraft at the NASA Goddard Space Flight Center and used for 
operational support of the Space Technology 5 (ST5) mission. 11-13 SpinKF2, introduced in Ref. 13, has been used to 
support the Time History of Events and Macroscale Interaction during Substorms (THEMIS) mission. 14 This paper 
introduces two new implementations, SpinKF3 and SpinKF4. In all the versions except SpinKFl, three of the 
components of the state error vector are the usual infinitesimal attitude error angles, making the usual 3x3 attitude 
covariance matrix a submatrix of the 6x6 covariance and greatly facilitating the computation of measurement 
sensitivity matrices. 
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The paper first presents the seven-parameter angular-momentum-based representation. This is followed by a brief 
introduction to SpinKF2 and complete derivations of SpinKF3 and SpinKF4. Data simulated to represent the 
THEMIS spacecraft are used to compare the performance of the four SpinKF variants and a more conventional EKF 
based on the quaternion and body rotation rate. Attitude estimation results using THEMIS flight data are presented 
and conclusions are drawn from the results of the simulations and flight data results. 


II. Angular-Momentum-Based Attitude Parameterization 

The spacecraft’s angular momentum about its center of mass expressed in its body frame, , or expressed in an 
inertial reference frame, L 7 , obeys the equations of motion 15 

dL B /dt = S B -( 0 BI x-L B (la) 

and 

dL I /dt = N I =A T BI K B (lb) 

as well as the constraint 

( 2 ) 

where N is the external torque, A m is the inertial-to-body attitude matrix, and the angular velocity (0 B1 is given by 




(3) 


In this equation J is the spacecraft moment of inertia tensor and h. nt is the angular momentum of any moving parts 

(reaction wheels, steerable antennas or solar arrays, flexible modes, fuel slosh, etc.) relative to the spacecraft. 
Spacecraft dynamics are commonly modeled by Eq. (la) and kinematics by 


dA BJ j dt - [© B1 x]A bi , 

where 


(4) 


[vx] = 


0 — v. 


-v 0 


3 V 2 
0 -v, 

v, 0 


(5) 


denotes the cross product matrix for an arbitrary 3 -component vector v. A quaternion or some other lower- 
dimensional representation of A m is often integrated rather than Eq. (4), but this distinction is not important for this 
paper. An alternative formulation uses Eq. (lb) and 


^ B ~ ^BI^I 


( 6 ) 


in place of Eq. (la). Both of these formulations have the disadvantage for application to spinning spacecraft that 
many components of the state vector are rapidly varying parameters, requiring small integration steps. 

The formulation in this paper is based on the observation that both the attitude matrix A BJ and the rotation matrix 


R B i = (L fl ■ L, )/ 3x3 - L,I4 + L*L' + (1 + L* • L, r l (L* * L, XL* x L, )' 


(7a) 


take the unit vector L 7 = L J L to = L B / L . Thus the product R r BI A BI takes L ; into itself, so it is a rotation 

about L j and can be expressed as 




(7b) 
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where 


P(e, 0) = (cos 0)/ 3x3 + (1 - cos 0)ee r - sin 0 [ex] (7c) 

is the matrix representing a rotation by angle 0 about axis e. Multiplying both sides of Eq. (7b) by R m gives 

Aft — Rfii R(Lj , £) . (8a) 

Multiplying this equation on the right by / 3x3 = P^ BI R BI and using the identity R BJ R(L I ,C ) )R r BI = R( L 5 ,£) gives 

A m =R(L B £)R Bn (8b) 

Either Eq. (8a) or (8b) gives our parameterization of the attitude matrix A m (\) in terms of the seven-component 
state vector 

X = [l/ B L Tj {J. (9) 

The dynamics are given by Eqs. (la), (lb), and a differential equation for £ derived in Ref. 1 and as Eq. (21) below. 
SpinKFl used a slightly different state vector, as described in the Appendix. 


III. Reduction to a Six-Component Error State 


A straightforward Kalman filter implementation would use the 7x7 covariance P = P{(Ax)(Ax) r } of the error 
vector Ax = x - x , where x = P{x} is the expectation of x. Note that an overbar will always denote an expectation. 
If the constraint of Eq. (2) were linear, P x would be singular. However, since Eq. (2) is a nonlinear constraint 
analogous to the unity norm constraint of the four-component quaternion, we expect P x to be ill-conditioned but not 

singular, in parallel with the situation observed in quaternion estimation. 9 As pointed out in Ref. 10, this means that 
all the useful covariance information can be contained in a well-conditioned matrix of rank six, in the same way that 
the quaternion norm constraint allows us to use a reduced-dimension covariance matrix. 2 In other words, in the 
context of extended Kalman filtering, we can treat P. as if it were singular, even though it is in fact only ill- 

conditioned. In the present case, the 7x7 covariance matrix P is assumed to have a null vector \ nu u that obeys 


^x X null ~ ^7x1 ' 


( 10 ) 


This is not really a null vector of the true covariance; it is a vector along the direction of the eigenvalue that leads to 
the ill-conditioning by being much smaller than the others. 

The null vector must be orthogonal to all the vectors representing physically possible errors in x. The six physically 
possible errors are a variation in £ two independent variations of L 5 perpendicular to h B , two variations of L 7 

perpendicular to L 7 , and simultaneous length-changing equal-magnitude variations of L, B in the direction of 
h B and of L ; in the direction of L 7 . Thus the normalized null eigenvector of P x must be 


X null 





( 11 ) 


which is a vector of errors violating the norm constraint of Eq. (2). 

We seek a 7x6 matrix S with a Moore-Penrose pseudoinverse 16 * S + such that 


and 


S + S = L 


6x6 5 


SS - ^7x7 X null X 


$ X null ~ ^6x1 * 


T 

null 9 


(12a) 

(12b) 

(12c) 
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Note that S and S + must be functions of expectations rather than true values, which are unknown by the estimator. 
The six-dimensional error vector 

Ay = *S + Ax (13) 

has a well-conditioned 6x6 error covariance P y given by 

P y * £{(Ay)(Ay) r } = S + P x (S + ) r . (14) 

Because of Eqs. (10) and (12b), the 7x7 error covariance P x can be recovered by 

SP y S T =SS + P x (S + ) T S T = SS + P X (SS + ) T = P x (15) 

and the seven-component error vector, which must be orthogonal to x null , by 

A* = 7 7X7 AX = ( ss+ + x nu// x L;;) Ax = SA y ■ ( 16 > 

Beginning with SpinKF2, the matrix S + has been chosen so that the first three components of Ay are equal to the 
three attitude error angles in the body reference frame, A0 . This choice results in the upper left-hand corner of P y 
being the usual 3x3 attitude error covariance, and we will see later that it also greatly facilitates computation of 
measurement sensitivity matrices. 

The expression for A0 as a function of Ax is found by computing the first-order increment A A m of A B/ (x) with 
variations AL 5 , AL ; , and A£ , and using the relationship 

[A0x] = -(AA BI )Agj . (17) 

This gives, after considerable algebra, 

A0= -ZT 1 t B x [AL b - A B1 (x) AL 7 ] + [ A£ - 1" 1 w • (AL^ + AL 7 )]L 5 , (18) 

where 

w=(l + L fl -L / r 1 (L s xL / ) = (I 2 + L b • L, )"' (L b x L ; ) . (19) 

Interpreting Eq. (18) as a relation of time variations of 0 , L B , L 7 , and £ gives 




dL, dC i 

+ i- L " 


After substituting Eq. (1), Eq. (20) reduces to the dynamic equation for £, 

d^jdt = (1 + t B ■ L J r 1 t(L B + 1, ) • (0 BI + L-\t B x L / ) • (N b + N 7 )] . (2 1) 

The upper three rows of S + are given by Eq. (18). The different varieties of SpinKF are distinguished by different 
choices for the lower three rows. 

A. SpinKF2 

SpinKF2 satisfies Eqs. (12) by choosing the lower three rows of S + to give 
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so Eq. (13) gives 


Ay 2 = 


A0 

al b + 'M*) AL /_ 


(23) 


Note our convention of using a numerical subscript n to indicate that an expression is valid only for a specific 
SpinKFfl algorithm. The appearance of A B/ (x) in Ay made SpinKF2 more complex than SpinKF3 and SpinKF4, 
so we will not present its further development here. The interested reader can find the details in Ref. 13. 

B. SpinKF3 


SpinKF3 chooses the lower three rows of S + to give 

> 1 {-[L B x]-t fl w 7 '} r , {[L fl x]/f B/ (x)-t B w r } 


s 3 + = 


T _1t it 
7 3x3 2 L B L 5 


It 


This means that the six-component error state vector is 

A0 


A y 3 


AL B - t(L b • AL b - L 7 • AL ; )L £ 


3x1 



' A0 1 


1 

£ 

to 

1 


(24) 


(25) 


because h B • AL B = AL = L 7 • AL ; for state errors that obey the angular momentum norm constraint. This Ay is the 
same as would be obtained in an EKF for a four-component quaternion and L 5 by using the techniques in Ref. 2 to 
handle the quaternion norm constraint. The pseudoinverse of S+ is 


s 3 = 


®3x3 

-4(x)[L b x] 


■■3x3 

A B,W 


l [£ b +t B x /< fl/ (x)w] r rw^oowf 


C. SpinKF4 

SpinKF4 chooses the lower three rows of S + to give 


5 4 + = 


t ~ If rr 

'B 

ir t t 

2 L, I L ‘B 


L 1 {-[L„x]-L„w 7 ’} L *{[L„x]^ g/ (x)-L fl w'} L 


T 


J B 1J i ^ * BI y 

1 — — r r T 

*3x3 2 L / L / 


0 


3x1 


and therefore 


A y 4 


A0 

AL I +Hi B -AL B -i I -AL I )t I 



' A0 " 


AL ; 


(26) 


(27) 


(28) 


The motivation for SpinKF4 is that AL ; varies more slowly than AL 5 , which should make the covariance 
propagation easier. The pseudoinverse of S+ is 


*4 = 


[L b x] 


^3x3 


A B i( x ) 


3x3 


(L 5 -L 5 xwf L V+^/Wwf 


(29) 
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IV. Kalman Filter Formulation 


A Kalman filter for the seven-component state vector x uses Eqs. (1) and (21) to propagate the state estimate 
between observations. The filter update for a measurement z = h(x) is given by 17 

x(+) = x(-) + K x [z- h(— )] , (30) 

where the arguments (-) and (+) denote estimates before and after the update, respectively, z denotes the measured 
value, h(-) = h(x(-)) , and K x is the Kalman gain. The gain is given by 

K x = P x (-)H t x [H x P x (-)H t x +R]-\ (31) 

where the measurement sensitivity matrix is 

H x = 3h(x)/3x (32) 

and R is the measurement error covariance. The covariance is updated by 

^(+) = ( / 7x7 -K X H X )P X (~), (33) 

To avoid using an ill-conditioned covariance matrix, we substitute Eq. (15) into Eq. (31), giving 

K x = S(-)K y , (34) 

where the matrix S must be evaluated with the pre-update estimate and K y is given by 

K y = Py (-)H T y [H y P y (-)H T y + R]~ l (35) 

with 

H y = H x S(r ) = [5h(x)/8x](9x/3y) = dh/dy . (36) 

A. Covariance Update 

Substituting Eq. (15) into Eq. (33) and using Eqs. (12a) and (12b) gives 

^+>=('«x«-VW->- (37) 

Because S(~) is evaluated with the pre-update estimates, both P y (~) and P y (+) are defined in the six-dimensional 
subspace of the seven-dimensional error state space orthogonal to x„ u //(-). We would like P y (+) to be defined in the 
six-dimensional subspace orthogonal to x ww// (+). This can be accomplished by mapping the covariance into the 

seven-dimensional space using S(~) and then down to the six-dimensional subspace using S + (+) , i.e. 

P y (++) = [S + (+)S(-)]P y (+)[S + (+)S(-)] T . (38) 


The matrix S + (+)S(-) is not orthogonal except in the limit of zero measurement update, where Eq. (12a) shows 

that it is equal to the identity. The transformation in Eq (38) is subtle, and it has been ignored in the past, 2 but 
Reynolds 18 has forcibly argued for its theoretical and practical importance in a similar context. Note that Eq. (38) 

marks the only appearance of the matrix S + in the final formulation of the EKF. The matrix S appears there and also 
in the state update. 

B. State Update: SpinKF3 

The state update in SpinKF3 is given by 


x(+) = x(-) + S 3 (-)^[z- h(-)] ^ x(-) + S 3 (-) 


A0 “ 
AL*/ 


( 39 ) 
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Specifically, 


Lfl(+) = L b (-) + AL fi 

L 7 (+) = £;(-)+ 4 (-)[AL b + A0 x L fl (-)] = Agj (-)[L B (-) + AL S + A0 X L a (-)] . 


(40a) 

(40b) 


The update preserves the norm constraint of Eq. (2) to first order but violates it in second order, in parallel with the 
quaternion case. 2 The cross-product term in Eq. (40b) appears to be the first-order approximation to a rotation, so we 
trust the magnitude of L B (+) more than that of L 7 (+) and renormalize the latter by 


L j (++) = |l b (+)| L j (+) 


(41) 


to restore the norm constraint. 

C. State Update: SpinKF4 

The state update in SpinKF4 is given by 


x(+) = x(-) + SA-)K\z - h(-)] = x(-) + SA-) 


A0 

AL, 


Specifically, 


L a (+) = L a H - A0 x L fl (— ) + A bi (-)AL 7 = A bi (-)[L 7 (-) + AL 7 ] - A0 x L a (-) 

l 7 (+) = l 7 (-)+al 7 , 


(42) 

(43a) 

(43b) 


which also violates the norm constraint in second order. Since the cross-product term occurs in Eq. (43a), the norm 
constraint in SpinKF4 is restored by 

L b (++) = |L 7 (+)|£ b (+). (44) 


D. Covariance Propagation 

The matrices F y and G y in the covariance propagation equation 

dP /dt = F P + P F T +G OG t 

y/ y y y y y^ y 

are most easily computed directly from the equations for the reduced state vector y. We assume that 


and 


N B = K B +Jn u 


h int= h m t - Jn v’ 


(45) 

(46a) 

(46b) 


where n w and n y are assumed to be independent zero-mean Gaussian white noise processes with 


E • 

n v (0 

\iri 

r 

■ = Q(t)8(t-t') = 

'Q v (0 

®3x3 


_ n „« 


j 


. ®3x3 

Q«(0_ 


(47) 


where 8(t-t') is the Dirac delta function. The definitions in Eq. (46) are arranged so that n w and n v have 
dimensions of angular acceleration and angular velocity, respectively, and Q u and Q v have units of rad 2 /sec 3 and 
rad 2 /sec. Note that Q v may be non-zero even for a nominally rigid spacecraft. 

To first order, the attitude error vector obeys the dynamics equation 19 

d( A0 )/ dt = A(O b1 - Sb B1 x A0 (48) 
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For SpinKF3 we need the dynamics equation obeyed by the angular momentum error in the body frame 
J(AL b )/ dt = A[N S - (0 B1 X L fl ] = L b x Am BI - m BI x AL B + Jn u . 

From Eq. (3) we have 

A(O m = J-\AL b - AL int ) = J-'AL B + n v , 

so differentiating Eq. (25) gives 


——=F y ,Ay, +Gy3 


with 


and 


F * m 


-[(B s/ X] J- 1 

. °3x3 [L B xy-'-[w B/ x]j 


G y 3 = 


^3x3 ®3x3 

[L„x] J 


SpinKF4 uses the dynamics equation obeyed by the angular momentum error in the inertial frame 
d(AL,)/dt = ANj = A[4(x)N b ] = A{4 (x)(/ 3)<3 + [A0 x])(N B + Jn„)} 


= ^bi ( X )(^ n „ - N B x A0) . 


We also use Eqs. (3) and (17) to obtain 


A<o b , = */ _1 A[^ s/ (x)L / — L /n( ] = J _i {-[A0x]^ b/ (x)L / + A B1 (x)AL, - AL jm } 


(49) 

(50) 

(51) 
(52a) 
(52b) 

(53) 


= J- 1 [L B xA6+^ / (x)AL / ] + n v . (54) 

An equation analogous to Eq. (51) is obtained for SpinKF4 with 

F = y'[L B x]-[M B/ x] J~'a bi (x) 

~ ^ b ;( x )[^ 5 X ] ® 3 x 3 

and 

G = hx3 °3x3 

_ |_®3x3 4^.' 

V. Filter Implementation 

SpinKF has been implemented in MATLAB™ as a subsystem of the Multimission Spin-Axis Stabilized Spacecraft 
Attitude Ground Support System that has supported NASA Goddard Space Flight Center missions for many years. 
The new EKF subsystem adds the capability to solve for a time-dependent attitude history and could be used for 
real-time applications, if needed. The software processes sensor data and presents it to the EKF as vector 
observations. After discarding outliers, the EKF integrates the state vector and its covariance to the next observation 
time using a 4 th -order Runge-Kutta integrator with an appropriate time step. The EKF obtains the spacecraft 
ephemeris and geomagnetic field and computes torques due to gravity gradients and any residual constant spacecraft 
magnetization at each integration step. 15 Control torques, if known, can be included. Then the sensor residual and 
the sensitivity matrix are computed, and the state and covariance are updated using the symmetric Joseph form 17 for 
the covariance update in place of Eq. (37). Tests were performed both with and without the transformation of Eq. 
(38). The transformation slowed the computations and did not provide consistently improved performance, so it was 
omitted in the numerical computations reported in this paper. 


(55a) 

(55b) 
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A. Singularity Avoidance 

It is clear throughout the development of this filter that the spacecraft angular momentum is required to be nonzero, 
and that the algorithm is also singular when L s and L 7 are 180° apart. The software checks for the latter singular 
condition and redefines the inertial reference frame so that L 7 is always greater than a user-specified distance from 
-L b in the modified frame, transforming all reference vectors along with L 7 . Reconstructing the attitude referenced 
to the standard inertial frame is only a matter of keeping track of these reference frame rotations, which are all 
handled internally and are totally transparent to the user of the software. 

B. Measurement Models 

It is easier to calculate the measurement sensitivity matrices directly from the six-component error vector y than 
from the seven-component state vector x. For a vector measurement, z = \ B and 

h ( x )= ^B/( x ) v / = : ( / 3 x3-[A0X]H/(x) v / =v b +v b xA0 , (56) 

where 

v B = A bi (x)\, . (57) 

It follows from Eq. (25) and the rightmost part of Eq. (36) that 

"T” °3x,]- < 58 > 

In the case of quaternion measurements, such as those output from an autonomous star tracker, an error quaternion is 
defined by 

Aq = [Aq T Aq 4 f =q(-)®q- 0 l bs , (59) 

where q(-) is constructed from the pre-update state x(-) , and where the quaternion multiplication convention of 
Ref. 2 is used, so that A(q®q) = A(q)A(q') . Then the observation residual is given by 

z - h(-) = -2 sign(Ag 4 ) SU1 ^ Aq (60) 

I Aq| 


and the measurement sensitivity matrix is simply 

H q ;a ,emion =[I^0^]. (61) 

For spacecraft with gyros, the gyro outputs are used as measurements to update the state, rather than being used for 
state propagation in model replacement mode. 2 The gyro sensitivity matrix is the only measurement sensitivity 
matrix that is different between SpinKF3 and SpinKF4. The gyro measurement model is 

h(x) = f/J~‘(L fi - L jn( ) , (62) 

where the rows of U are the gyro sensitive axes in the body frame. Thus in SpinKF3 we have simply 

//J° = t/[0 3X 3 J~ 1 ]. (63) 

In SpinKF4, on the other hand, Hf ro = U[J~ X 0 3x3 0 3xl ] , which yields 

Hf™ - Hf ro S 4 (-) = UJ~ X [[L fl x] 4 b/ (x)]. (64) 
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VI. Testing with Simulated Data 


Several tests were performed to compare the performance of different versions of SpinKF and the Unit Vector 
Filter 20 (UVF) that has supported many three-axis stabilized spacecraft at Goddard Space Flight Center over the past 
16 years. For application to spinning spacecraft, the UVF was modified to estimate the quaternion and rotation rate 
rather than the gyro biases, to use dynamics propagation rather than gyro propagation, and to use linearized 
dynamics matrices for covariance propagation that are related to that of SpinKF3 by 

F 

1 UVF 

and 

-l 

G y 3 . (65b) 

A. Simulation Parameters 

The simulator allows modeling of spinning spacecraft with a wide variety of sensor types and several torque 
scenarios. The tests exercised all the key features of the filters using simulation parameters based on the THEMIS 
series of spinning spacecraft. 14 THEMIS consists of five identical probes designed to study magnetic substorms in 
the Earth’s magnetosphere. The probes will be maneuvered into orbits with periods and phasing that cause them to 
line up radially from the Earth at apogee once every four days. Apogee will occur in the magnetotail during the first 
part of the primary science campaign from December 2007 to April 2008. 

THEMIS was launched on February 17, 2007. The initial orbit for all five probes had an inclination of about 14°, 
period of 31 hours, and eccentricity of 0.85. The perigee and apogee heights were roughly 900 km and 87000 km, 
respectively. The spin axes have all since been oriented along the normal to the ecliptic plane, but the simulations 
use the initial release spin direction that had a right ascension of 5 1 .3° and declination of -23.2°. 

The probes are low mass (126 kg), spin- stabilized at 20 revolutions per minute (rpm) about the body Z-axis, and 
carry single-head slit-type Sun sensors and three-axis magnetometers (TAM) for attitude determination. The slit Sun 
sensor measures the Sun angle from the Z-axis once per spin and generates a timing pulse. The ground software 
constructs a body frame Sun vector from the measured angle and the known slit azimuth. The TAM measurement 
frequency used for attitude determination is 8 Hz, which gives a generous 24 observations per spin period. 

The TAM is both an attitude sensor and a science instrument, and its inherent noise is very small (0.01 nT). Thus, 
the TAM error for attitude determination is almost entirely due to errors in the reference field, which were modeled 
in the simulation and in the filters as Gaussian errors with standard deviation of 100 nT on each axis. The 
International Geomagnetic Reference Field (IGRF) was used to model the magnetic field. The Sun sensor errors are 
due to a combination of Sun angle measurement error and Sun pulse timing error. They were modeled as Gaussian 
errors with standard deviation of 0.16° per axis. 

The inertia tensor used for all the simulations is 


J= diag([13, 13,22]) kg-m 2 

(66) 

This value is of the correct magnitude for THEMIS after deployment of the TAM boom but prior to deployment of 
the long wire radial booms. The tensor was taken to be axisymmetric to simplify computation of simulated attitudes. 
A nutation angle of 2° was assumed in all the torque-free runs. The simulator and the filter propagation routine 

allow for environmental torques, but these were turned off for the tests presented here. 
The process noise was taken to be 


Q v = 10" 6 rad 2 /sec x diag([l, 1, 3]) 

(67a) 

and 


Q u = 1CT 7 rad 2 /sec 3 x diag([l, 1, 3]) 

(67b) 


GjjVF ~ 


hx3 ®3x3 
^3x3 J 


hx3 ®3x3 
®3x3 J 


"|-i 




hx3 ^3x3 
_^3x3 J _ 


(65a) 
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The initial state covariance matrix was constructed assuming an initial attitude uncertainty of 20° about X and Y and 
1 80° about Z and a rate uncertainty of 10 deg/sec on each axis. 

The simulation time span is 25 minutes and is centered on perigee where the altitude is low enough for the TAM 
data to be useful. The first three minutes of each run are discarded from the results to allow time for the filter 
transients to damp out. 

B. Test Results 

The values reported in this section are a measure of the spin axis pointing error, which is the root-sum-square of the 
X- and Y - axis errors. The spin phase error is not included since this is much less important operationally. The Z- and 
T-errors are obtained by converting the filter state vector estimates into a history of equivalent quaternions, 
comparing with the truth model quaternions, and determining the root-mean-square error about each axis. 

A series of tests were designed to stress the filters in four different ways: reducing the number of observations per 
spin period, increasing the propagation step size, increasing the initial attitude error, and increasing the error in the 
inertia tensor model. 

Prior to these stressing cases, four baseline scenarios were tested to verify filter performance with parameters similar 
to those used for operational support of THEMIS. The Nominal scenario consists of torque-free motion, a 2° 
nutation angle, no sensor misalignments, an accurately modeled inertia tensor, and uses TAM data (24 per spin 
period) and Sun sensor data (one per spin period). The second scenario was the Nominal scenario with gyro data 
added. The third scenario was the same as the Nominal scenario, but with misalignments of 0.1° on the TAM and 
Sun sensor. The fourth scenario was also similar to the Nominal, but with a torque applied in the body frame for 
10% of each spin period to precess the spin axis by 30°. Note that the filter is NOT given the torque values for any 
of the tests presented here; if the torque values are passed to the filter, the filter will agree very closely with the truth 
model (the errors then are indistinguishable from the nominal case). 

Table 1 shows the results of these first tests for SpinKFl, 2, 3, 4, and the modified UVF described by Eq. (65). In 
these tests, the initial attitude errors relative to the truth model are 10° about the X - and T-axes, 45° about Z, and the 
rate errors are 5 deg/sec about Zand Y, and 10 deg/sec about Z. 

It is clear from Table 1 that all five filters are performing well and are in good agreement. As a measure of 
computational burden, the mean clock times for these runs are 47, 29, 25, 25, and 26 sec for SpinKFl, SpinKF2, 
SpinKF3, SpinKF4, and the UVF, respectively. These times are platform-dependent, but the relative times are 
significant. SpinKFl is slower in part because its partial derivative expressions are much more complicated than the 
other filters, but also because its MATLAB™ code has not been as thoroughly optimized. It is not clear why 
SpinKFl performs better than the other filters with the attitude slew here and in the next example. 


Table 1. Results of baseline filter tests. 


Scenario 

SpinKFl 
Pointing Error 

SpinKF2 
Pointing Error 

SpinKF3 
Pointing Error 

SpinKF4 
Pointing Error 

UVF 

Pointing Error 

Nominal 

0.0810° 

0.0873° 

0.0873° 

0.0874° 

0.0873° 

Nominal plus gyro 

0.0769° 

0.0249° 

0.0249° 

0.0249° 

0.0249° 

Misaligned sensors 

0.1325° 

0.1388° 

0.1387° 

0.1390° 

0.1392° 

30° Slew 

(no torque data to KF) 

0.1142° 

0.7224° 

0.7139° 

0.7348° 

0.7110° 


The next set of tests examines how the filters respond with reduced data availability. Table 2 shows results for the 
same four scenarios as in Table 1, but with the TAM data frequency reduced from 24 to 2.1 observations per spin 
period. (The number is deliberately nonintegral so the spin phase at the TAM measurement times will not be 
constant.) The propagation step size was set to 1/8 second to match the step size in the baseline tests (driven in that 
case by the TAM frequency). 
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Table 2. Filter test results with TAM data frequency reduced to 2.1 observations per spin period. 


Scenario 

SpinKFl 
Pointing Error 

SpinKF2 
Pointing Error 

SpinKF3 
Pointing Error 

SpinKF4 
Pointing Error 

UVF 

Pointing Error 

Nominal 

0.1171° 

Diverged 

0.1097° 

0.1098° 

0.1097° 

Nominal plus gyro 

0.1120° 

0.0364° 

0.0367° 

0.0370° 

0.0367° 

Misaligned sensors 

0.1596° 

Diverged 

0.1538° 

0.1537° 

0.1538° 

30° Slew 

(no torque data to KF) 

0.1625° 

1 .2726° 

1 .2662° 

1 .2654° 

1.2764° 


With fewer data, the attitude errors shown in Table 2 are larger than those in Table 1, but with the exception of 
SpinKF2, the filters still perform well. When the TAM frequency is reduced further, all the filters tend to diverge, 
depending largely on the chance transients from the first few data points. 

The filter determines the propagation step size first, from the data time steps, and second, from an input parameter 
that specifies a maximum allowed time step. Whenever the time between observations is larger than this maximum 
step size, additional integration steps are inserted as required. To test the numerical propagation, the next test 
increased the maximum allowed step size while keeping the actual data frequency low. With the TAM frequency at 
2.1 observations per spin period and the maximum step size equal to 0.5 sec, all the filters except SpinKFl and 
SpinKF4 begin to fail the test scenarios. Table 3 shows that increasing the step size further to 0.7 sec causes 
SpinKF4 to fail in some cases, but displays no clear pattern. The convergence region is small for all the filters when 
the data are sparse and the propagation step size is large. Whether a filter diverges or not then is largely just a 
happenstance caused by the first few data points. If the transients push the state far enough from the truth so that 
nonlinearities in the sensor models become important, the filter often never recovers. When the step size is increased 
to 1 sec, SpinKFl also diverges. 


Table 3. Filter test results with TAM data frequency reduced to 2.1 observations per spin period and the 

maximum propagation step size increased to 0.7 sec. 


Scenario 

SpinKFl 
Pointing Error 

SpinKF2 
Pointing Error 

SpinKF3 
Pointing Error 

SpinKF4 
Pointing Error 

UVF 

Pointing Error 

Nominal 

0.7027° 

Diverged 

Diverged 

Diverged 

Diverged 

Misaligned sensors 

0.7115° 

Diverged 

Diverged 

0.1816° 

Diverged 

30° Slew 

(no torque data to KF) 

0.1631° 

Diverged 

1 .4762° 

1 .2690° 

Diverged 


The next series of tests considers the effect of varying the initial attitude error in the nominal scenario (torque-free, 
2° nutation, 8 Hz TAM, no misalignments, 1/8 sec propagation step size). The initial rate error is kept at 5 deg/sec 
about X and Y, and 10 deg/sec about Z. Table 4 shows selected results. Entries with relatively large values indicate 
cases where the filter did converge but took longer than three minutes to settle down to linear behavior, so some of 
the initial transients are included in the statistics. The results are mixed, but it is clear the SpinKF filters do not 
outperform the UVF when presented with large initial errors. 

The last row in Table 4 shows that all five filters converged when the X and Y errors were zero and the Z error was 
large, that is, with spin phase error only. The initial spin phase typically is completely unknown. Experience with 
many runs with simulations and actual flight data shows that occasional filter divergence can be expected even when 
the initial pointing direction is known to within a few degrees. When this happens, restarting the filter with the phase 
shifted by 180° almost always gives a good solution. To avoid these occasional unnecessary divergences, the filter 
could self-initialize by making a rough attitude estimate using data from the first few spin periods. This feature will 
be added to the operational version of the filter in a future release. 
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Table 4. Filter results for nominal scenario with various initial attitude errors. 


Initial Attitude Errors 

SpinKFl 

Pointing 

Error 

SpinKF2 

Pointing 

Error 

SpinKF3 

Pointing 

Error 

SpinKF4 

Pointing 

Error 

UVF 

Pointing 

Error 

Xerr = 30°, Yerr = 30°, Zerr = 0 

0.0802° 

0.0873° 

0.0873° 

0.0874° 

0.0873° 

Xerr = 60°, Yerr = 60°, Zerr = 0 

0.0802° 

0.0873° 

3.1932° 

0.1119° 

0.0873° 

Xerr =80°, Yerr = 80°, Zerr = 0 

Diverged 

0.0873° 

14.4281° 

Diverged 

0.1199° 

Xerr = 90°, Yerr = 0, Zerr = 0 

0.0808° 

0.0873° 

0.0873° 

0.0874° 

0.0873° 

Xerr = -90°, Yerr = 0, Zerr = 0 

0.0810° 

16.9711° 

0.0891° 

0.0874° 

0.0873° 

Xerr = 1 20°, Yerr = 0, Zerr = 0 

0.0810° 

2.2748° 

Diverged 

0.0875° 

0.0874° 

Xerr = 0, Yerr = 1 20°, Zerr = 0 

Diverged 

Diverged 

Diverged 

Diverged 

Diverged 

Xerr = 1 70°, Yerr = 0, Zerr = 0 

11.2017° 

Diverged 

0.0873° 

Diverged 

Diverged 

Xerr = -170°, Yerr = 0, Zerr = 0 

3.7612° 

Diverged 

Diverged 

0.0876° 

Diverged 

Xerr = 30°, Yerr = 30°, Zerr = -175° 

0.0804° 

16.9142° 

0.0833° 

Diverged 

0.0815° 

Xerr = 0, Yerr=0, Zerr =175° 

0.0807° 

0.0890° 

0.0875° 

0.0881° 

0.0873° 


The final set of tests is again based on the nominal scenario, but now includes various amounts of error in the inertia 
tensor. The simulator uses the inertia given by Eq. (66) in all cases, but the filter here is given an incorrect inertia for 
propagation. The first row in Table 5 uses an inertia tensor close to the actual value for THEMIS after TAM boom 
deployment. In this case, all the filters except SpinKFl show similar errors of only 1.1°. The errors are relatively 
small here because the inertia tensor error consists primarily of a rotation of the principal axis frame about the 
nominal body frame Z-axis by about 41° with only 0.15° X and Y rotation. The 41° Z rotation has no significant 
effect since the truth model inertia tensor is axisymmetric. The combined X and Y rotation tips the major principal 
axis away from the body Z-axis, which induces coning motion into the propagation that is not present in the truth 
model attitude. However, a larger source of error in this example is that this tensor is far from axially symmetric, so 
the nutational motion is not predicted accurately. The eigenvalues of J are approximately 7.8, 18.5, and 21.7 kg-m 2 ; 
the lack of symmetry is indicated by the difference between the first and second eigenvalues not being small relative 
to the third. 

In the second case in Table 5, the principal axis frame is tipped by 6.6° from the body frame. Compared with the 
first case, the symmetry error is less important since the eigenvalues are 11.0, 15.3, and 21.8 kg-m 2 . Thus, the 
coning error predominates. In the third and fourth cases in Table 5, the principal axis frame is tipped by 10.3° and 
14.1°, respectively, and again, coning error predominates. As the coning angle increases, the SpinKF4 error follows 
it closely, but the other filters begin to show larger errors. The covariance propagation is probably more accurate for 
SpinKF4 since it uses AL/ rather than AL 5 in its error state vector. Specifically, the lower half of F y 4 in Eq. (55a) 
does not include a factor of J, so this part of the covariance propagation is not corrupted by the error in the inertia 
tensor. 
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Table 5. Filter results for nominal scenario with various errors in the inertia tensor. 


Inertia Tensor 
(kg-m 2 ) 

SpinKFI 

Pointing 

Error 

SpinKF2 

Pointing 

Error 

SpinKF3 

Pointing 

Error 

SpinKF4 

Pointing 

Error 

UVF 

Pointing 

Error 

^ EKF ~ 

" 12.4 - 5.3 -. 0 : 
- 5.3 14.0 -.0 
-.02 -.01 21.7 


1 .9872° 

1.1143° 

1.1134° 

1.1164° 

1.1137° 

7 ekf — 

1 


7.4997° 

6.5419° 

7.0813° 

6.8028° 

7.3736° 

7 ekf ~ 

' 14 0 - 1 . 5 " 

0 12 0 
- 1.5 0 22 


11.9418° 

10.3834° 

14.5196° 

10.5923° 

13.7981° 

7 ekf ~ 


1 

15.9104° 

17.8481° 

30.1164° 

13.9457° 

38.6044° 


VII. Estimation Results with Flight Data 

The SpinKF2 version of the filter has been used for ground attitude support for the five THEMIS probes. Before 
filter performance could be assessed, the sensors needed to be calibrated. In particular, the TAM was deployed after 
launch on a 2-meter boom and has an alignment uncertainty of approximately 1°. The calibration utility uses an 
attitude-independent method 21 to determine biases, scale factors, and skewness of the axes, and an iterative method 22 
to determine the orthogonal sensor alignment. SpinKF2 is a key part of the iterative method. The calibrations were 
repeated several times for each probe using data sets from March, 2007. Each data span covers 50 to 60 minutes 
near perigee. The reference magnetic field was obtained using the IGRF 2005 model to 10 th order along with 
definitive ephemerides for the spacecraft positions. 

The spacecraft attitude history was determined using SpinKF2 before and after calibration. Before calibration, the 
sensor residuals from the filter show offsets of 1° to 2° from zero primarily due to misalignment. Calibration 
removes these offsets and reduces their standard deviations by up to a factor of two. 

The spin vector is nominally aligned with the body Z-axis and is constant in an inertial frame; however, a nonzero 
coning angle is expected, especially prior to deployment of the radial wire booms. The calibration reduces the 
standard deviation of the Z-axis pointing direction by roughly a factor of two, yielding a post-calibration 3-sigma Z- 
axis uncertainty of 0.4° to 0.5° for Probes A-D and 0.7° for Probe E. This is in excellent agreement with the Z-axis 
coning motion predicted from the inertia tensors. The angular momentum direction is more nearly constant; it does 
not vary due to coning. The angular momentum direction 3-sigma uncertainties range from less than 0.1° to 0.3° for 
all post-calibration solutions for the five probes. 

Figure 1 shows a typical solution for the estimated angular momentum right ascension, declination, and magnitude 
for Probe A. The variation in the angular momentum direction seen in Figure 1 is not physical. Gravity gradient is 
the largest environmental torque, and it precesses the angular momentum vector by only 0.01° per orbit, much less 
than the variation seen in Figure 1 . It is likely that this small variation is due to error in the reference field caused by 
a small, but not entirely negligible, ephemeris error. This problem was seen, greatly exaggerated, early in the 
mission before accurate ephemerides were available. On the other hand, the apparent slow downward drift of the 
angular momentum magnitude results from cooling of onboard fuel that, being forced to the outer edge of the tanks, 
contracts away from the spin axis, thereby increasing the moment of inertia about that axis and decreasing the spin 
rate. Assuming a constant inertia tensor, the filter misinterprets this as a decreasing angular momentum. 
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Fig. 1. Right ascension, declination, and magnitude of the angular momentum vector estimated using 
SpinKF2 for THEMIS Probe A after calibration of the magnetometer. 

VIII. Conclusions 

A class of extended Kalman filters for spinning spacecraft has been developed using a seven-component angular- 
momentum-based attitude parameterization. The state vector is constrained by the requirement that the magnitude of 
the angular momentum is the same in all reference frames. This constraint allows the filters to employ a six- 
component error state instead of the error vector of the full seven-component state, in parallel with the procedure 
commonly used to estimate the constrained four-component quaternion representation of attitude. Different filters in 
this class all use the same seven-component state but differ in their specification of the six-component error state. 
Two of these filters have been presented previously, and two new formulations are introduced in this paper. The six- 
component error states of the two filters introduced here comprise the vector of infinitesimal attitude error angles 
and the angular momentum vector in either the spacecraft body frame or in an inertial frame. 

These four filters, along with a conventional quaternion-based filter, were tested using data simulated to represent a 
spacecraft of the THEMIS mission. A nominal test was performed for torque- free motion, with a two degree 
nutation angle, no sensor misalignments, an accurately modeled inertia tensor, three-axis magnetometer data, and 
Sun sensor data. Three variations on this nominal test were run by adding gyro data, misaligning the magnetometer 
and Sun sensor by 0.1°, and applying a torque to precess the spin axis by 30°. The filter was not given the torque 
values for any of the tests because that would have resulted in estimation errors indistinguishable from the nominal 
case. The performance of the filters was virtually identical in this series of tests, with the exception that the oldest 
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version of the filter, SpinKFl, performed significantly better in the 30° slew case. The reason for this is unknown, 
and the execution time of this filter was almost twice as great as that of the two filters introduced in this paper, in 
part because its partial derivative expressions are much more complicated, but also because its code has not been as 
thoroughly optimized. The five filters would have all performed equally well in the 30° slew case if they had been 
provided with perfect torque information. 

This baseline set of cases was followed by a series of tests designed to stress the filters in four different ways: 
reducing the number of observations per spin period, increasing the propagation step size, increasing the initial 
attitude error, and increasing the error in the modeled moment-of-inertia tensor. The oldest version of the filter, 
SpinKFl, and the new variant, SpinKF4, employing the infinitesimal attitude error angles and the angular 
momentum components in an inertial reference frame as the error state showed the most robust performance in these 
tests. On balance, the SpinKF4 version of the filter showed the best combination of robustness and efficiency in the 
simulations. 


Appendix: SpinKFl 

SpinKFl uses the seven- component state vector 

X - [L; L T b Cj. (Al) 

This differs from Eq. (9) in using the unit vector L 7 in place of L ; , and in an unimportant transposition of 
components. The constraint on this state vector is not Eq, (2), but |l ; | = 1 . Thus the null vector is 

*„„//= ±[4 0i x4 J. (A2) 


Equations (12)— (13) are satisfied by 


1 

0 

0 

K, 

®2x4 

0 

1 

0 


®4x3 


^4x4 . 


(A3) 


where 


^=( z . L / )/ 3x 3 


-L ; z r +z\J I +(l + z-L / ) _I (zxL / )(zxL / ) r 


(A4) 


is a rotation matrix that takes L y into z = [0 0 1 ] 7 . The six-dimensional EKF follows the same pattern as the other 
versions of SpinKF, except that the state renormalization after the measurement update is equivalent to 


L/(++) = L / (+)/|l / (+) . 


(A5) 


The computation of the measurement sensitivity matrices in SpinKFl is more difficult than in the other versions of 
SpinKF because the infinitesimal attitude error angles are not part of the reduced error vector. SpinKFl also 
contains a singularity when L/ and z are 180° apart that is not present in the other filters. This singularity is handled 
in the same way by redefining the inertial reference frame. SpinKFl is described in detail in Ref. 10. 
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